Interval Truth Model

Visualizations of Priors in the ITM

Authors
Affiliation

Matthias Kloft

University of Marburg

Björn S. Siepe

University of Marburg

Daniel W. Heck

University of Marburg

Published

July 21, 2025

1 Latent Appraisal

1.1 True Intervals: Informative Beta Prior on Marginal Locations and Widths

1.1.1 Sample From Priors

Show the code
# priors on marginal locations and widths
df <- data.frame(Tr_wid_splx = rbeta(1e5, 1.2, 3))
# multiplicative shifting parameters
s <- rbeta(1e5,1,1)
# lower bounds
df$Tr_L_splx <- 
  (1 - df$Tr_wid_splx) * s
# transform from simplex to bivariate normal
Tr_bvn <- simplex_to_bvn(
  cbind(
    df$Tr_L_splx,
    df$Tr_wid_splx,
    df$Tr_L_splx + df$Tr_wid_splx))

df <- cbind(df, Tr_loc = Tr_bvn[, 1], Tr_wid = Tr_bvn[, 2])

# transform back to simplex
Tr_splx <- bvn_to_simplex(
  cbind(df$Tr_loc, df$Tr_wid))
df$Tr_loc_splx <- Tr_splx[,1] + .5 * Tr_splx[,2]
df$Tr_wid_splx <- Tr_splx[,2]

1.1.2 Marginal Distributions

Show the code
df %>%
  ggplot() +
  geom_histogram(
    aes(Tr_loc_splx),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Tr_wid_splx),
               fill = "purple",
               alpha = .3,
               binwidth = .005) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

Show the code
df %>%
  ggplot() +
  geom_histogram(
    aes(Tr_loc),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Tr_wid),
               fill = "purple",
               alpha = .3,
               binwidth = .005, ) +
  scale_x_continuous(
    limits = c(NA, NA),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

1.1.3 Joint Distribution

1.1.4 Ternary

Show the code
df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Tr_loc_splx, y = Tr_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

1.1.4.1 Unbounded Scale

Show the code
df %>%
  ggplot() +
  geom_point(
    aes(x = Tr_loc, y = Tr_wid),
    size = 1,
    alpha = .1,
    shape = 16
  ) +
  theme_pubr()

1.2 Lambda (Item Discernibility)

1.2.1 sigma_lambda

Show the code
df$sigma_lambda <- exp(rnorm(1e5, log(.5), .5))
#median(lambda$sigma_lambda)

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_lambda),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.2.2 lambda

Show the code
df$lambda_loc <- exp(rnorm(1e5, 0, df$sigma_lambda))
df$lambda_wid <- exp(rnorm(1e5, 0, df$sigma_lambda))

df %>%
  ggplot() +
  geom_histogram(aes(x = lambda_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3 E (Person Proficiency)

1.3.1 mu_E

Then narrow Student t prior ensures that the population proficiency is regularized to 1.

Show the code
df$mu_E <- rnorm(1e5)

df %>%
  ggplot() +
  geom_histogram(aes(x = exp(mu_E)),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.2 sigma_E

Show the code
df$sigma_E <- exp(rnorm(1e5, log(.5), .5))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_E),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.3 E

Show the code
df$E_loc <- exp(rnorm(1e5, mean = df$mu_E, sd = df$sigma_E))
df$E_wid <- exp(rnorm(1e5, mean = df$mu_E, sd = df$sigma_E))

df %>%
  ggplot() +
  geom_histogram(aes(x = E_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.3.4 Lambda and E joint

Show the code
df %>%
  ggplot() +
  geom_histogram(aes(x = 1 / lambda_loc / E_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "SD of Latent Appraisal", y = "Density") +
  theme_pubr()

1.4 Omega (Residual Correlation)

Show the code
beta_param <- 2
df$omega <- rbeta(1e5, shape1 = beta_param, shape2 = beta_param)*2-1

df %>%
  ggplot() +
  geom_histogram(aes(x = omega),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(-1, 1), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

1.5 Latent Appraisals

1.5.1 Marginal Appraisals

Show the code
A <- rbvnorm(1e5, 
             mean1 = df$Tr_loc, 
             mean2 = df$Tr_wid,
             sd1 = 1 / df$lambda_loc / df$E_loc,
             sd2 = 1 / df$lambda_wid / df$E_wid,
             cor = df$omega)

df$A_loc <- rnorm(1e5, mean = df$Tr_loc, sd = 1 / df$lambda_loc / df$E_loc)
df$A_wid <- rnorm(1e5, mean = df$Tr_wid, sd = 1 / df$lambda_wid / df$E_wid)

A_splx <- bvn_to_simplex(df[,c("A_loc", "A_wid")])
df$A_loc_splx <- A_splx[,1] + .5 * A_splx[,2]
df$A_wid_splx <- A_splx[,2]

df %>%
  ggplot() +
  geom_density(aes(A_loc_splx),
               fill = "lightblue",
               alpha = .3,
               shape = 16) +
  geom_density(aes(A_wid_splx),
               fill = "purple",
               alpha = .3,
               shape = 16) +
  scale_x_continuous(limits = c(0, 1), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Appraisal Locations / Widths", y = "Density") +
  theme_pubr()

1.6 Joint Appraisals

Show the code
df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = A_loc_splx, y = A_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

2 Biases

2.1 Shifting Bias

2.1.1 sigma_b

Show the code
df$sigma_b <- exp(rnorm(1e5, log(.5), 1))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_b),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

2.1.2 b

Show the code
df$b_loc <- rnorm(1e5, mean = 0, sd = df$sigma_b)
df$b_wid <- rnorm(1e5, mean = 0, sd = df$sigma_b)

df %>%
  ggplot() +
  geom_histogram(aes(x = b_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(-5, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

Show the code
b_splx <- bvn_to_simplex(df[,c("b_loc", "b_wid")])
df$b_loc_splx <- b_splx[,1] + .5 * b_splx[,2]
df$b_wid_splx <- b_splx[,2]

df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = b_loc_splx, y = b_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

2.2 a (Scaling Bias)

2.2.1 sigma_a

Show the code
df$sigma_a <- exp(rnorm(1e5, log(.5), .5))

df %>%
  ggplot() +
  geom_histogram(aes(x = sigma_a),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 2.5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

2.2.2 a

Show the code
df$a_loc <- exp(rnorm(1e5, mean = 0, sd = df$sigma_a))

df %>%
  ggplot() +
  geom_histogram(aes(x = a_loc),
                 fill = "lightblue",
                 binwidth = .01) +
  scale_x_continuous(limits = c(0, 5), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  theme_pubr()

3 Responses

3.1 Marginal Responses

Show the code
df$Y_loc <- df$A_loc * df$a_loc + df$b_loc
df$Y_wid <- df$A_wid + df$b_wid

Y_splx <- bvn_to_simplex(df[,c("Y_loc", "Y_wid")])
df$Y_loc_splx <- Y_splx[,1] + .5 * Y_splx[,2]
df$Y_wid_splx <- Y_splx[,2]

df %>%
  ggplot() +
  geom_histogram(
    aes(Y_loc_splx),
    fill = "lightblue",
    alpha = 1,
    binwidth = .005
  ) +
  geom_histogram(aes(Y_wid_splx),
               fill = "purple",
               alpha = .3,
               binwidth = .005) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Marginal Locations / Widths", y = "Density") +
  theme_pubr()

3.2 Joint Responses

Show the code
df %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Y_loc_splx, y = Y_wid_splx),
    size = 1,
    alpha = .1,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()


4 Recovery Check

4.1 Define Parameter Bounds for Variances

Show the code
# compute a benchmark for the mean and SD of the parameters
mean_benchmark <- simplex_to_bvn(c(.4, .2, .4), type = "ilr")
sd_benchmark_loc <- simplex_to_bvn(c(.98, .01, .01), type = "ilr")
sd_benchmark_wid <- simplex_to_bvn(c(.495, .01, .495), type = "ilr")
# mean for Tr_loc
mu_Tr_loc <- mean_benchmark[1]
# mean for Tr_wid
mu_Tr_wid <- mean_benchmark[2]
# SD forTr_loc
sigma_Tr_loc <- sd_benchmark_loc[1] / 4
# SD Tr_wid
sigma_Tr_wid <- abs(sd_benchmark_wid[2] - mean_benchmark[2]) / 4

# SDs for other parameters
sigma_lambda_E_loc <- .3
sigma_lambda_E_wid <- .3
sigma_a_loc <- .3
sigma_b_loc <- sigma_Tr_loc / 3
sigma_b_wid <- sigma_Tr_wid / 3

4.2 Generate Data for One Replication

Show the code
n_respondents <- 50
n_items <- 20

sim_data <-
  generate_itm_data_sim_study(
    n_respondents = n_respondents,
    n_items = n_items,
    mu_Tr_loc = mu_Tr_loc,
    mu_Tr_wid = mu_Tr_wid,
    sigma_Tr_loc = sigma_Tr_loc,
    sigma_Tr_wid = sigma_Tr_wid,
    sigma_lambda_E_loc = sigma_lambda_E_loc,
    sigma_lambda_E_wid = sigma_lambda_E_wid,
    sigma_a_loc = sigma_a_loc,
    sigma_b_loc = sigma_b_loc,
    sigma_b_wid = sigma_b_wid,
    omega_beta = 3
  )
responses <- sim_data$responses

4.2.1 Plot Data

Show the code
samples <- gather_values(
  lower = responses$x_L,
  upper = responses$x_U,
  item_id = responses$jj,
  step_size = 0.01
)
samples %>% 
  ggplot(aes(x = samples)) +
  geom_density(fill = "purple", alpha = .5) +
  facet_wrap(~item_id) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

4.2.2 Plot True Intervals as Error Bars

Show the code
data.frame(
  idx = 1:n_items,
  lower = sim_data$parameters$Tr_L,
  upper = sim_data$parameters$Tr_U
) %>%
  ggplot(aes(y = idx)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper), height = 0.3) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
  labs(x = "Latent Truth",
       y = "Item") +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

4.2.3 Plot True Intervals as Ternary Plot

Show the code
data.frame(
  idx = 1:n_items,
  Tr_L = sim_data$parameters$Tr_L,
  Tr_U = sim_data$parameters$Tr_U) %>% 
  mutate(
  Tr_loc = Tr_L + (Tr_U - Tr_L) / 2,
  Tr_wid = Tr_U - Tr_L
) %>%
  ggplot() +
  geom_abline(intercept = 0, slope = 2) +
  geom_abline(intercept = 2, slope = -2) +
  geom_point(
    aes(x = Tr_loc, y = Tr_wid),
    size = 4,
    alpha = .5,
    shape = 16
  )+
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
    scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  theme_pubr()

4.3 Fit Model

4.3.1 Stan Data

Show the code
### Stan data declaration
I <- length(unique(responses$ii))
J <- length(unique(responses$jj))
N <- nrow(responses)
ii <- responses$ii
jj <- responses$jj
nn <- c(1:N)
Y_splx <- cbind(responses$x_splx_1, responses$x_splx_2, responses$x_splx_3) %>%
  as.matrix()

#Y_splx[is.na(Y_splx)]

# sanity check: all rows sum to 1
# check <- apply(Y_splx, 1, sum) %>% as.data.frame()
# table(check)

## stan data list
stan_data <- list(
  I = I,
  J = J,
  N = N,
  ii = ii,
  jj = jj,
  nn = nn,
  Y_splx = Y_splx,
  link = 1
)

4.3.2 Compile Model

Show the code
# Choose model to fit
model_name <- "itm_simulation_v2_beta"
# Compile model
model <-
  cmdstanr::cmdstan_model(
    stan_file = here("src", "models", paste0(model_name, ".stan")),
    pedantic = TRUE,
    quiet = FALSE
  )

4.3.3 Run Sampler

Show the code
# number of MCMC chains
n_chains <- 4
# Run sampler
fit_prior_check <- model$sample(
  data = stan_data,
  seed = 2023,
  chains = n_chains,
  parallel_chains = n_chains,
  iter_warmup = 500,
  iter_sampling = 500,
  refresh = 500,
  adapt_delta = .99,
  init = .1
)
Running MCMC with 4 parallel chains...

Chain 1 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 2 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 3 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration:   1 / 1000 [  0%]  (Warmup) 
Chain 4 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 4 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 1 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 1 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 2 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 2 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 3 Iteration: 500 / 1000 [ 50%]  (Warmup) 
Chain 3 Iteration: 501 / 1000 [ 50%]  (Sampling) 
Chain 4 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 4 finished in 213.4 seconds.
Chain 2 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 2 finished in 218.9 seconds.
Chain 1 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 1 finished in 222.3 seconds.
Chain 3 Iteration: 1000 / 1000 [100%]  (Sampling) 
Chain 3 finished in 263.2 seconds.

All 4 chains finished successfully.
Mean chain execution time: 229.5 seconds.
Total execution time: 263.4 seconds.
Show the code
# save fit_prior_check
fit_prior_check$save_object(
  file =  here("fits", paste0(model_name, "_fit_prior_check.rds")))
Show the code
# load  fit
fit_prior_check <- readRDS(
  file =  here("fits", paste0(model_name, "_fit_prior_check.rds")))

4.3.4 Get Estimates

Show the code
estimates_summary <- fit_prior_check$summary()

4.3.5 Sampler Diagnostics

Show the code
fit_prior_check$diagnostic_summary()
$num_divergent
[1] 0 0 0 0

$num_max_treedepth
[1] 0 0 0 0

$ebfmi
[1] 0.6690134 0.7145780 0.6690593 0.6188122
Show the code
# sampler diagnostics
fit_prior_check$sampler_diagnostics(format = "df") %>% 
  psych::describe(quant = c(.05,.95),) %>%
  round(2) %>%  
  as.data.frame() %>% 
  dplyr::select(median, min, Q0.05, Q0.95,  max) %>% 
  .[-c(7:9),]
              median    min  Q0.05  Q0.95    max
treedepth__     7.00   6.00   6.00   7.00   7.00
divergent__     0.00   0.00   0.00   0.00   0.00
energy__      540.30 469.52 502.22 579.02 624.73
accept_stat__   0.99   0.91   0.97   1.00   1.00
stepsize__      0.04   0.04   0.04   0.05   0.05
n_leapfrog__  127.00  63.00 127.00 127.00 255.00
Show the code
# convergence diagnostics
convergence_summary <- 
  fit_prior_check$draws(format = "df") %>%
  summarise_draws(.x = ., "rhat", "ess_bulk", "ess_tail") %>%
  remove_missing() %>%
  select(-variable) %>%
  psych::describe(., quant = c(.05, .95)) %>%
  as.data.frame() %>%
  select(mean, median, Q0.05, Q0.95, min, max)

convergence_summary %>% round(3)
             mean   median    Q0.05    Q0.95     min      max
rhat        1.002    1.002    1.000    1.005   0.999    1.010
ess_bulk 1956.513 1861.096  938.414 3327.645 417.098 4357.585
ess_tail 1467.217 1476.239 1161.061 1756.521 682.151 2000.306

4.3.6 Effective sample size (ESS) & Rhat Plots

Show the code
# color scheme
color_scheme_set(scheme = "purple")

# Effective sample sizes
plot_neff <-
  mcmc_neff_hist(bayesplot::neff_ratio(fit_prior_check), binwidth = .01) +
  labs(title = "A") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  )
# Rhat
plot_rhat <-
  bayesplot::mcmc_rhat_hist(bayesplot::rhat(fit_prior_check)) +
  labs(title = "B") +
  guides(color = FALSE, fill = FALSE) +
  theme(
    legend.text = element_blank(),
    legend.key = element_blank(),
    title = element_text(size = 16, face = "bold")
  ) +
  yaxis_text(on = TRUE)

# Combined plot
plot_diagnostics <- gridExtra::grid.arrange(plot_neff, plot_rhat, ncol = 2)

4.4 Parameter Plots

Show the code
# color scheme
color_scheme_set(scheme = "purple")

true_parameters <- sim_data$parameters

4.4.1 Person Competence: Location

Show the code
plot_E_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("E_loc"), true = true_parameters$E_loc) +
  labs(
    subtitle = "Person Competence: Location",
  ) 
plot_E_loc

4.4.2 Person Competence: Width

Show the code
plot_E_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("E_wid"), true = true_parameters$E_wid) +
  labs(subtitle = "Person Competence: Width")
plot_E_wid

4.4.3 Person Scaling Bias: Location

Show the code
plot_a_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("a_loc"), true = true_parameters$a_loc) +
  labs(subtitle = "Person Scaling Bias: Location") +
  theme(
    text = element_text(
      family = "serif",
      size = 16,
      colour = 1
    ),
    plot.title = element_text(size = 20,
                              face = "bold"),
    plot.margin = unit(c(.1, .1, .1, .1), "cm"),
    axis.title = element_text(size = 16, color = 1),
    axis.title.x = element_text(margin = margin(t = 5)),
    axis.title.y = element_text(margin = margin(r = 5)),
    axis.text = element_text(size = 12, colour = "black"),
    axis.text.x = element_text(margin = margin(t = 3))
  )
plot_a_loc

4.4.4 Person Shifting Bias: Location

Show the code
plot_b_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("b_loc"), true = true_parameters$b_loc) +
  labs(subtitle = "Person Shifting Bias: Location") 
plot_b_loc

4.4.5 Person Shifting Bias: Width

Show the code
plot_b_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("b_wid"), true = true_parameters$b_wid) +
  labs(subtitle = "Person Shifting Bias: Width")
plot_b_wid

4.4.6 Item Discernability: Location

Show the code
plot_lambda_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("lambda_loc"),
                       true = true_parameters$lambda_loc)
plot_lambda_loc

4.4.7 Item Discernability: Width

Show the code
# plot
plot_lambda_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("lambda_wid"),
                       true = true_parameters$lambda_wid) +
  labs(subtitle = "Item Discernability: Width")
plot_lambda_wid

4.4.8 Residual Correlation

Show the code
plot_omega <-
  mcmc_recover_scatter(x = fit_prior_check$draws("omega"),
                       true = true_parameters$omega) +
  xlim(-1, 1) +
  ylim(-1, 1) +
  labs(subtitle = "Residual Correlation")
plot_omega

4.4.9 Latent Truth Intervals

Show the code
# plot
plot_Tr_loc <-
  mcmc_recover_scatter(x = fit_prior_check$draws("Tr_loc"),
                       true = true_parameters$Tr_loc) +
  labs(subtitle = "Item True Location")
plot_Tr_loc

Show the code
# plot
plot_Tr_wid <-
  mcmc_recover_scatter(x = fit_prior_check$draws("Tr_wid"),
                       true = true_parameters$Tr_wid) +
  labs(subtitle = "Item True Width")

plot_Tr_wid

Show the code
# ITM estimates
latent_truth_est_ilr <- data.frame(
  idx = 1:J,
  type = "estimated",
  Tr_L_est = fit_prior_check$summary("Tr_splx") %>% 
    as.data.frame() %>% 
    pull(median) %>% 
    .[1:J],
  Tr_wid_est = fit_prior_check$summary("Tr_splx") %>% 
    as.data.frame() %>% 
    pull(median) %>% 
    .[(J + 1):(J * 2)]
) %>%
  mutate(Tr_loc_est = Tr_L_est + Tr_wid_est / 2, 
         Tr_U_est = Tr_L_est + Tr_wid_est)

# compute means in the unbounded space
means_ilr <-
  cbind(jj, simplex_to_bvn(Y_splx)) %>%
  as.data.frame() %>%
  dplyr::group_by(jj) %>%
  dplyr::summarise(across(everything(),mean, na.rm = TRUE)) %>%
  dplyr::ungroup() %>%
  dplyr::select(-jj) %>%
  bvn_to_simplex() %>%
  mutate(
    idx = 1:J,
    type = "simple_mean",
    Tr_L_sm = x_1,
    Tr_U_sm = 1 - x_3
  ) %>%
  select(-c(x_1, x_2, x_3))


# true parameters
latent_truth_true <- data.frame(
  idx = 1:J,
  type = "true",
  Tr_L_true = true_parameters$Tr_L,
  Tr_U_true = true_parameters$Tr_U
)

df_interval_plot_ilr <- full_join(latent_truth_est_ilr, latent_truth_true) %>%
  full_join(means_ilr) %>% 
  mutate(type = factor(type, levels = c("estimated", "true")))

4.4.9.1 Interval Plot

Show the code
cols <- c("True" = "grey70","ITM (ILR)" = "red", "Mean" = "blue")

df_interval_plot_ilr %>%
  ggplot(aes(y = idx)) +
    geom_errorbarh(
    aes(xmin = Tr_L_true, xmax = Tr_U_true, col = "True"), 
    height = 0,
    linewidth = 5
  ) +
  geom_errorbarh(
    aes(xmin = Tr_L_sm, xmax = Tr_U_sm, col = "Mean"),
    height = 0,
    linewidth = 3,
    alpha = .4
  ) +
  geom_errorbarh(
    aes(xmin = Tr_L_est, xmax = Tr_U_est, col = "ITM (ILR)"),
    height = 0,
    linewidth = 2,
    alpha = .5
    ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = seq(0, 1, .25),
    expand = expansion()
  ) +
    scale_y_discrete(
    expand = expansion(add = 1)
  ) +
  scale_color_manual(values = cols) +
  labs(x = "Latent Truth", y = "Item") +
  theme_pubr() +
  theme(plot.margin = margin(.1, .5, .1, .1, "cm"),
        panel.grid.major = element_line())

5 Prior vs. Posterior

5.1 True Intervals

Show the code
latent_truth <- data.frame(
  Tr_loc_splx_post = fit_prior_check$draws("Tr_loc_splx") %>%
    unlist() %>%
    as.vector(),
  Tr_wid_splx_post = fit_prior_check$draws("Tr_wid_splx") %>%
    unlist() %>%
    as.vector()
) %>%
  mutate(
    jj = factor(rep(1:J, each = nrow(.) / J)),
    Tr_wid_splx_prior = rbeta(nrow(.), 1.2, 3),
    Tr_loc_splx_prior = (1 - Tr_wid_splx_prior) * rbeta(nrow(.), 1, 1) + Tr_wid_splx_prior / 2
  )

Tr_bvn <- simplex_to_bvn(
  cbind(
    latent_truth$Tr_loc_splx_prior - .5 * latent_truth$Tr_wid_splx_prior,
    latent_truth$Tr_wid_splx_prior,
    1 - (
      latent_truth$Tr_loc_splx_prior + .5 * latent_truth$Tr_wid_splx_prior
    )
  )
)
latent_truth <- cbind(latent_truth,
                      Tr_loc_prior = Tr_bvn[, 1], Tr_wid_prior = Tr_bvn[, 2])


latent_truth %>%
  ggplot() +
  geom_abline(intercept = 0,
              slope = 2,
              alpha = .7) +
  geom_abline(intercept = 2,
              slope = -2,
              alpha = .7) +
  geom_density_2d_filled(
    aes(x = Tr_loc_splx_prior, 
        y = Tr_wid_splx_prior),
    binwidth = .05, 
    alpha = .3, ) +
  geom_point(
    aes(x = Tr_loc_splx_post, y = Tr_wid_splx_post, col = jj),
    size = .3,
    alpha = .3,
    shape = 16
  ) +
  scale_x_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    labels = c("0", ".25", ".5", ".75", "1"),
    expand = expansion()
  ) +
  labs(x = "Location", y = "Width") +
  # supress legends
  guides(col = FALSE, fill = FALSE) +
  theme_pubr()

5.2 Hyper-Priors

5.2.1 sigma_lambda: Location

Show the code
colors <- c("prior" = "lightblue", "posterior" = "purple")

sigma_lambda_loc <- data.frame(posterior = fit_prior_check$draws("sigma_J[1]", format = "list") %>% unlist() %>%  as.vector())
sigma_lambda_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_loc %>%
  ggplot() +
  geom_density(aes(prior), fill = "lightblue", alpha = 1) +
  geom_density(aes(exp(posterior * .5 + log(.5))), fill = "purple", alpha = .3) +
  scale_x_continuous(limits = c(NA, 3), expand = expansion()) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.2 sigma_lambda: Width

Show the code
sigma_lambda_wid <- data.frame(
  posterior = fit_prior_check$draws("sigma_J[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_lambda_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_lambda_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.3 sigma_E: Location

Show the code
sigma_E_loc <- data.frame(
  posterior = fit_prior_check$draws("sigma_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_loc$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.4 sigma_E: Width

Show the code
sigma_E_wid <- data.frame(
  posterior = fit_prior_check$draws("sigma_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
sigma_E_wid$prior <- exp(rnorm(nrow(sigma_lambda_loc), log(.5), .5))

sigma_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(exp(posterior*.5+log(.5))),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(NA, 3),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.5 mu_E: Location

Show the code
mu_E_loc <- data.frame(
  posterior = fit_prior_check$draws("mu_I[1]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_loc$prior <- rnorm(nrow(sigma_lambda_loc))

mu_E_loc %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(posterior),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()

5.2.6 mu_E: Width

Show the code
mu_E_wid <- data.frame(
  posterior = fit_prior_check$draws("mu_I[2]", format = "list") %>% unlist() %>%  as.vector()
)
mu_E_wid$prior <-  rnorm(nrow(sigma_lambda_loc))

mu_E_wid %>% 
  ggplot() +
  geom_density(
    aes(prior),
    fill = "lightblue",
    alpha = 1  ) +
  geom_density(aes(posterior),
               fill = "purple",
               alpha = .3
               ) +
  scale_x_continuous(
    limits = c(-4, 4),
    expand = expansion()
  ) +
  scale_y_continuous(limits = c(0, NA), expand = expansion()) +
  labs(x = "Prior / Posterior", y = "Density") +
  scale_color_manual(values = colors) +
  theme_pubr()